{-# LANGUAGE TypeFamilies #-}
{- |
Module      :  Interpolation
Description :  <optional short text displayed on contents page>
Copyright   :  (c) Sergey Sichevskiy 2012
License     :  BSD3

Maintainer  :  s.sichevskij@gmal.com
Stability   :  experimental
Portability :  portable

-}


module Interpolation ( 
   lerp,
   biLerp,
   triLerp,
   lerpR,
   linearSampling,
   linearSamplingU,
 ) where

import Helper             ( floor' )
import GHC.Float          ( int2Double )

import Data.AdditiveGroup ( (^+^), (^-^) )
import Data.VectorSpace   ( VectorSpace, Scalar, (*^), )
import Data.Array         ( listArray, (!) )
import Data.Array.Unboxed ( UArray )

import qualified Data.Array.Unboxed as U

-- |Функция 'order' упорядочивает ...
lerp' :: (VectorSpace v, Floating a, RealFrac a,  a ~ Scalar v) => (Int -> v) -> (a -> v)
--lerp' g x =  f `par` (h `pseq` (f ^+^ k *^ h))
lerp' g x | abs k <= 1e-10 = f
          | otherwise      = f ^+^ k *^ (f' ^-^ f)
  where 
--    h  = f `par` (f' `pseq` (f' ^-^ f))
    f  = g $ x'
    f' = g $ x' + 1
    x' = floor x
    k  = x - fromIntegral x'

-- |Linear interpolation
lerp :: (VectorSpace v, Double ~ Scalar v) => (Int -> v) -> (Double -> v)
--lerp g x = f `par` (f' `pseq` (f ^+^ k *^ h))
lerp g x | abs k <= 1e-10 = f
         | otherwise      = f ^+^ k *^ (f' ^-^ f)
  where
--    h  = f `par` (f' `pseq` (f' ^-^ f))
    f  = g $ x'
    f' = g $ x' + 1
    x' = floor' x
    k  = x - int2Double x'

-- |Bilinear interpolation
biLerp :: (VectorSpace v, Scalar v ~ Double) => (Int -> Int -> v) -> Double -> Double -> v
biLerp f  = lerp . lerp f

-- |Trilinear interpolation
triLerp :: (VectorSpace v, Scalar v ~ Double) => (Int -> Int -> Int -> v) -> Double -> Double -> Double -> v
triLerp f  = biLerp . lerp f

class Interpolation a where
  type R a :: *
  lerpR :: a -> R a


instance ( Int ~ a
         , Interpolation b , VectorSpace b
         , Double ~ (Scalar b)
         ) 
      =>  Interpolation ( a -> b ) where
   type R (a -> b) = Double -> R b
   lerpR g = lerpR.lerp g


instance Interpolation Double where
   type R Double = Double
   lerpR a = a

-- |
linearSampling :: (Enum a, Ord a, RealFrac a, Floating a, Num b, VectorSpace b, a ~ Scalar b) => (a,a) -> a -> [(a,b)] -> (a -> b)
linearSampling (l,r) h xs = lerp' f . (\x -> (x-l)/h)
  where
    f i | i < 0 || i > n = 0
        | otherwise = arr ! i

    n   = floor $ (r-l)/h 
  
    arr = listArray (0,n) $ map snd $ sampling (l,r) h xs

linearSamplingU :: (Double,Double) -> Double -> [(Double,Double)] -> (Double -> Double)
linearSamplingU (l,r) h xs = lerp f . (\x -> (x-l)/h)
  where
    f i | i < 0 || i > n = 0
        | otherwise = arr U.! i

    n   = floor' $ (r-l)/h 

    arr :: UArray Int Double
--    arr = U.listArray (0,n) $ smoothGauss 10 80 $ map snd $sampling (l,r) h xs
    arr = U.listArray (0,n) $ map snd $ sampling (l,r) h xs

--sampling :: (Double,Double) -> Double -> [(Double,Double)] -> [(Double,Double)]
sampling :: (Enum a, Ord a, Fractional a, Num b, VectorSpace b, a ~ Scalar b) => (a,a) -> a -> [(a,b)] -> [(a,b)]
sampling (l,r) h vs = loop [l,l+h .. r] $ zip vs $ tail vs
 where
  loop [] _ = []
  loop (x:xs) [] = (x,0):loop xs []
  loop (x:xs) v@((v0,v1):vs) | x < x0 = (x,0):loop xs v
                             | x > x1 = loop (x:xs) vs
                             | otherwise = (x,t) : loop xs v
    where (x0,y0) = v0
          (x1,y1) = v1
--          t = y0 ^+^ (y1^-^y0) ^* (x-x0)/(x1-x0)
          t = y0 ^+^ k *^ (y1^-^y0)
          k = (x-x0) / (x1-x0)


{--
smoothGauss :: Double -> Int -> [Double] -> [Double]
smoothGauss sigma krad xs =  U.elems arr
  where
    arr :: UArray Int Double
    arr = U.accumArray (+) 0 (0,n) $ filter (\(i,_) -> i >= 0 && i <= n) $ concat $ zipWith kernel xs [0..]

    kernel v i = [ (x + i, v * gauss x) | x <- [-krad .. krad] ]

    gauss x = 1 / (sigma * sqrt( 2*pi )) * exp( -0.5 * (fromIntegral x / sigma)**2 )

    n = length xs - 1
--}